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Abstract 

The Widom-Rowhnson model of a fluid mixture is studied using a new cluster 
algorithm that is a generalization of the invaded cluster algorithm previously 
applied to Potts models. Our estimate of the critical exponents for the two- 
component fluid are consistent with the Ising universality class in two and 
three dimensions. We also present results for the three-component fluid. 
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Some years ago Widom and Rowlinson introduced a simple but non-trivial continuum 
model that exhibits a phase transition 0]. The two-component formulation of this model 
consists of "black" and "white" particles; particles of the same type do not interact, but 
particles of differing type experience a hard-core repulsion at separations less than or equal 
to cr 0. In this letter we present a new Monte Carlo method for simulating the Widom- 
Rowlinson (WR) model and apply the method to study the properties of the transition point 
in two and three dimensions. 

There have been relatively few Monte Carlo studies of the WR critical point because 
of the combined difficulties of treating hard-core systems and critical slowing down using 
standard Monte Carlo techniques. The algorithm presented here overcomes these difficulties 
by using a generalization of the cluster approach introduced by Swendsen and Wang [Q. 
The algorithm is a variant of the invaded cluster (IC) method [§,0 adapted to continuum 
problems. We find that the algorithm has almost no critical slowing down and that we can 
obtain accurate values of the critical density and the exponent ratios (3/v and 7/z/ with 
relatively modest computational effort. 

The 2-component WR model is expected to be in the Ising universality class. Our results 
for f3/v and 7/1/ are consistent with this assumption and are in good agreement with the best 
known values for the Ising model. Our value for the critical density of the three-dimensional 
{d = 3) WR model agrees with recent results obtained by Shew and Yethiraj [0. 

We also consider a WR model in which there are q components, any two of which interact 
via a hard-core repulsion. Despite the apparent similarity to the g-state Potts model, the 
phase structure may have additional features Our algorithm easily extends to these 

g-component WR models, and we present results for the 3-component model in d = 2, 3. 

Graphical Representations of the WR Model. A configuration of the WR fluid consists 
of two sets of points, S and T, corresponding to the positions of the black and white 
particles respectively. In the grand canonical ensemble, the probability density for finding 
the configuration {S, T) is 
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Z is the grand partition function, Zx (^2) is the fugacity of the black (white) particles, and A^i 
(A^2) is the number of black (white) particles. The object F, which expresses the hard-core 
interaction between particles of different types, vanishes if any point in S is within a distance 
a of any point in T and is one otherwise. Symmetry considerations insure that the critical 
point is along the line z^ = Z2; hereafter we restrict attention to the case, z = zi = Z2. 

To motivate and justify our cluster algorithm, we consider a different representation of 



the WR model. In the "gray" representation IT^jO one considers the particles without 
reference to color, but with configurations weighted according to the number of allowed 
colorings. Let be a list of points, W = (ri, . . . rjy). Clusters of particles can be defined 
by the condition that every particle in a cluster is within a distance a of some other particle 
in the cluster. Particles in a cluster must all be the same color, so if there are C{W) distinct 
clusters (including single particles), there are 2^^^^ allowed colorings. Starting with Eq. (|T]) 
and working through the combinatorics, we find that the probability density for W is given 
by 

PiW) = ^^2^^^). (2) 

These densities describe the gray measures. The appropriate gray measure for the q- 
component model is defined by the analog of Eq. (|]) with 2 replaced by q. To return 
to the distribution in Eq. (|l|) starting from the gray representation, we select one of the 
2C{W) ^Qj, qC(w)-^ allowed colorings with equal probability. It turns out that the gray mea- 
sures are a special case of the models studied in ||12[. The connection with the WR models 
was discussed in |T3[ and made precise in [0,|11[ . 



Cluster Algorithms. With the gray representation in mind, we consider the following 
cluster algorithm. Starting from a configuration W of gray particles, clusters are identified. 
Each cluster is independently labeled black or white with probability 1/2 and all the white 
particles are removed. In the next step, white particles are replaced via a Poisson process 



at fugacity z in the free volume permitted by the black particles. In the final step, color 
identities are erased and we obtain the next gray configuration. The generalization to the 
g-component WR model is that the fraction of clusters deleted is This algorithm is the 
generalization of the Swendsen-Wang method to the WR model. The algorithm is described 
in more detail and detailed balance is proved in and independently in [|T^]. 



Because we are interested in efficiently sampling the transition point of WR models, we 
will forsake a Swendsen-Wang approach in favor of an IC algorithm. The steps involving the 
coloring and discarding of clusters are essentially the same, but rather than repopulating 
the free volume by a fixed fugacity process, particles are sequentially added with a uniform 
distribution throughout the free volume until a stopping condition is fulfilled. For example, 
one could add particles until a fixed particle number is reached. If a stopping rule is chosen 
that enforces a condition which is characteristic of criticality, then a critical state of the 
system is sampled automatically. In a finite volume the IC method samples an ensemble 
that differs from the canonical ensemble, but it presumably converges to the correct infinite 
volume distribution for all local observables. The validity of the IC method and its relation 
to the Swendsen-Wang method is discussed briefly below and in detail in Ref. 

The signature of the phase transition in the WR model is percolation of a gray clus- 



ter [0,|r^. Thus an appropriate stopping rule for the IC algorithm is the spanning of a 
gray cluster - particles are added at random in the allowed volume until a cluster spans the 
system. (In our case, we use periodic boundary conditions and spanning is said to occur 
when a cluster wraps the torus.) The other modiflcation is that the spanning cluster is 
erased on each deletion move, thereby ensuring that a new spanning cluster can form in the 
repopulation move. 

If A^tot denotes the total number of points that are needed to satisfy the stopping con- 
dition, then z = Niot/V is an estimator for the critical fugacity Zc- In the limit V —>■ oo, 
it is reasonable to assume that the distribution for z becomes sharp. If this assumption is 
valid, each move is identical to a move of the Swendsen- Wang-type algorithm at the peak 
value of z. It follows that this peak value is Zc - no other value of the fugacity would exhibit 



a critical cluster. Hence, if the distribution for z is very narrow in a finite volume, the IC 
algorithm is essentially the Swendsen- Wang-type algorithm with small fluctuations in the 
fugacity. 

Results. We now present our simulation results using the IC method. We collected 
statistics for the following quantities: the average number of particles in the spanning cluster, 
M; the normalized autocorrelation function TM{t) of the number of particles in the spanning 
cluster as a function of "time" t as measured in Monte Carlo steps; the compressibility x 
defined by 

X=^T.sl (3) 

where Sj is the mass of the ith cluster and the spanning cluster is included in this sum; 
the estimator for the critical fugacity z = (A^tot)/^; the fluctuations in z, a1 = ((A^tot) ~ 
{Ntot)'^)/V'^; and the average number of gray particles per unit area (volume) p which is an 
estimator of the critical density. System size is measured in units of the particle diameter, 
a. The results for these quantities are presented in Tables 1-lV along with our estimate of 
the statistical errors. 

The critical exponents that depend on the magnetic exponent t/h can be obtained from 
the fractal dimension of the spanning cluster via M ~ or from the compressibility via 
X ~ L'^^'^. The various critical exponents are related hj D = y^, 7/1^ = 2yh — d, and 
/3/z/ = d — yh- The dynamical properties of the algorithm can be measured by the integrated 
autocorrelation time defined by tm = \ + Ylu=i ^Mif)- The integrated autocorrelation time 
is roughly the number of Monte Carlo steps between statistically independent samples and 
enters into the error estimate for M. In practice, it is necessary to cut off the upper limit of 
the sum defining tm when Tm becomes comparable to its error. The increase in tm defines 
a dynamic exponent zm via tm ~ . 

Our results for the 2-component, d = 2 WR fluid are summarized in Table |. From the 
log-log plot of M versus L shown in Fig. we find that yh = d — [3 / u = 1.873 ± 0.002, and 
hence /3/z/ ~ 0.127, a value consistent with the exact Ising result of = 1/8. Similarly, a 
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log-log plot of X versus L yields ^jv = 1.743 ± 0.003, consistent with the exact Ising value, 
7/z/ = 7/4. These results support the hypothesis that the 2-component d = 2 WR fluid is 
in the Ising universality class. The errors quoted here are associated with the least squares 
fitting procedure and are two standard deviations. We also estimated the statistical errors 
in Uh and 7/z/ by generating synthetic data sets consistent with the estimated errors in the 
measured values of M and x ^"^^ found similar results. 

From the data of Table | we have estimated the infinite volume critical values of the 
density, pc and fugacity, Zc. A linear fit for p(L) versus 1/L yields pc = 1.5652. A three 
parameter fit of the form 

p{L) = - A/L^ (4) 

yields pc = 1.5662 with x = 0.96. Because the biggest source of error is the uncertainty 
in the fitting form rather than the statistical errors in the raw data, we estimate the error 
in Pc as several times the difference between these two fits. Hence, we conclude that pc = 
1.566 ±0.003. 

Within the statistical error the fugacity z is unchanged for the three largest system sizes. 
We take these values and several times the statistical error to estimate the critical fugacity, 
Zc = 1.726 ± 0.002. To our knowledge, there are no independent estimates of pc and Zc for 
the 2-component, d = 2 WR fluid. 

The autocorrelation function, TM{t) decreases rapidly and oscillates about zero after 
t ^ 10. Our results for tm for various system sizes are summarized in Table |. The slow 
increase of tm with L indicates that the dynamic critical exponent zm is small or zero 
{tm ~ InL). Because of its small value and our limited data, we cannot make a more 
precise statement. Fitting rj\,/(t) to a single exponential leads to decorrelation times similar 
to the integrated autocorrelation time. 

Our results for the 2-component, d = 3 WR fluid are summarized in Table 0. Power 
law flts of M and x versus L yield y,^ = 3~ f3/u = 2.479 ± 0.001 and 7/// = 1.961 ± 0.003, 
respectively. These values are consistent with each other and with the recent estimate of 
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i/h = 2.4815(15) obtained in Ref. |T6[. The results confirm the expectation that the 2- 
component, d = 3 WR model is in the d = 3 Ising universality class. 

In Fig. ^ we show p versus 1/L. A linear fit yields pc = 0.7484. A three parameter fit 
of the form given in Eq. yields pc = 0.7478 with x = 1.16. We estimate the error in pc 
as several times the difference between these two fits and conclude that pc = 0.748 ± 0.002. 
This value of pc is in agreement with and improves upon the recent result in Ref. 0, 
Pc = 0.762 ± 0.016. These values for pc are much higher than older estimates of pc which 
were in the range of 0.41 to 0.57 []T7 |. 



From the estimates of tm shown in Table |I|, we see that tm does not appear to increase 
with L. It may be that there is no critical slowing for IC dynamics for the 2-component, 
d = 3 WR model as is the case for the d = 3 Ising model under IC dynamics . 

The fluctuations in the estimator of the fugacity, a^, decrease with L as ~ with 
a ^ 0.5 for d = 2 and a ~ 0.8 for d = 3. The d = 2 value of a is the same as was found 
for the d = 2 Ising model in the IC ensemble while for = 3 it is somewhat larger than the 



results obtained for the d = 3 Ising model |TB[ where a = 0.69 ± 0.01. The fact that cr^ 



as the system size increases insures that the IC ensemble is close to the canonical ensemble. 

Our results for the 3-component WR fluid in d = 2 are summarized in Table |ITT[ Using 
the same finite size scaling analysis we used for the 2-component WR fluid, we find that 
D = Uh = 1.842 ± 0.004 and hence (3/v ~ 0.16. This result for Uh is consistent with 
our observed value of 7/z/ = 1.681 ± 0.008. These results are not consistent with the 
corresponding value for the 3-state, d = 2 Potts model where yh = 28/15. Even more 
surprising, our estimated value of yh for the WR fluid is less than the minimum value of yh 
for any d = 2 Potts model with a continuous transition {yh ~ 1.86603 for q = 3.332). This 
observation deserves further study. Is the 3-component, d = 2 WR fluid outside the Potts 
universality classes or are the systematic errors considerably larger than we have guessed? 

We also find that p approaches pc as 1/L and extract the value pc = 1.657 ± 0.001. The 
standard deviation of the estimator of the fugacity decreases as ~ A log-log plot 

of Tm versus L yields the estimate of the dynamic exponent zm = 0.58, that is, for this case 



z is sufficiently large for us to conclude that 2; > 0. 

Our results for the 3-component, d = 2 WR fluid are summarized in Table |IV|. The 
corresponding 3-state, = 3 Potts model is believed to have a ffist-order transition, and it is 
likely that this behavior holds for the 3-component d = 3 WR fluid. On the other hand, both 
the observed values of M and x ^"^^ well-described by power laws. This situation also holds 
for the 3-state Potts model for computationally accessible system sizes and reflects the fact 
that the transition is very weakly first-order. More study is needed to determine the order of 
the transition for this case of the WR model. The IC method finds the transition temperature 
for Potts models independently of whether the transition is ffist-order or continuous [§. 
Hence, we believe that extrapolated values of z yields the transition value of the fugacity, 
Zf. 1.16, and that p ~ 0.795 lies between the density of the two co-existing phases at the 
transition. The standard deviation of the fugacity decreases with L as (Tj; ~ The 
effective dynamic exponent describing the increase in tm is zm = 0.62. 

We have shown that cluster methods may be effectively used to study the Widom- 
Rowlinson model. To our knowledge, the IC algorithm is the ffist example of an algorithm 
that performs efficiently near the critical point of a continuum system with hard-core inter- 
actions. We have obtained accurate values of the critical density and fugacity for Widom- 
Rowlinson models in d = 2 and 3. The 2-component Widom-Rowlinson model appears to 
be in the Ising universality class. However, the 3-component d = 2 model deserves further 
study and might not be in the 3-state Potts universality class. 

This work was supported by NSF grants DMR-9632898 and DMR-9633385. 
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TABLES 

TABLE I. Dependence of M, p, x, z, az, and tm on L for the 2-component, d = 2 WR fluid. 

The error estimates represent one standard deviation. The averages are over 10^ spanning clusters. 



L M 




P 


X 


z 






40 1511(1) 




1.5247(4) 


1584(2) 


1.7201(7) 


0.212 


0.58 


60 3233(3) 




1.5379(3) 


3217(5) 


1.7247(6) 


0.170 


0.60 


80 5527(5) 




1.5450(3) 


5289(9) 


1.7267(6) 


0.150 


0.72 


120 11836(12) 




1.5516(3) 


10760(20) 


1.7265(5) 


0.120 


0.78 


160 20282(20) 




1.5552(2) 


17752(30) 


1.7262(4) 


0.101 


0.77 


TABLE IL Dependence of M, p, x, z, a^, 


and Tm on L for the 2-component, d = 


3 WR fluid. 


The averages are over 


10^ 


spanning clusters. 










L M 




P 


X 


z 






10 313.0(1) 




0.74022(7) 


119.5(2) 


0.9387(1) 


0.138 


0.59 


20 1745.4(6) 




0.74440(4) 


466.1(9) 


0.940(1) 


0.077 


0.57 


30 4768(2) 




0.74567(2) 


1031(2) 


0.9403(1) 


0.056 


0.57 


TABLE IIL Dependence of M, p, x, z, Uz 


, and Tm on L for the 3-component, d = 


2 WR fluid. 


The averages are over 


10^ 


spanning clusters. 










L M 




P 


X 


z 




TM 


40 1480(2) 




1.6124(8) 


1543(4) 


1.965(3) 


0.364 


0.88 


80 5325(7) 




1.6352(6) 


4987(13) 


1.960(1) 


0.280 


1.2 


120 11211(18) 




1.6418(6) 


9810(30) 


1.953(1) 


0.237 


1.7 


160 19013(32) 




1.6455(6) 


15869(50) 


1.949(1) 


0.213 


1.9 
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TABLE IV. Dependence of M, x, z, cr^, and tm on L for the 3-component, d = 3 WR fluid. 

The averages are over lO'' spauuing cluslers. 



L 


M 


P 


X 


z 




Tm 


10 


299.6(3) 


0.7914(2) 


109.8(2) 


1.1789(8) 


0.230 


0.58 


20 


1639(2) 


0.7947(2) 


409.1(9) 


1.1717(6) 


0.145 


0.74 


30 


4407(6) 


0.7947(1) 


874(2) 


1.1670(5) 


0.117 


0.83 
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FIGURES 

FIG. 1. Log-log plot of M, the average mass of the spanning cluster, versus L, the lin- 
ear dimension of tlic lattice for the 2-componcnt, d = 2 WR fluid. A least squares fit yields 
D = 2-l3/v = 1.873. 

FIG. 2. Plot of p versus 1/L for the 2-component, d = 3 WR fluid. A least squares fit yields 
Pc = 0.7484. 
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